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PREFACE 


This  report  is  intended  to  be  a  brief  summary  of  the  most 
basic  elements  of  the  sublet  of  Power  Spectral  Analysis  of  time- 
series  data.  These  elements  are  presented  and  discussed  heuristlcally 
without  rigorous  mathematical  justification.  It  is  hoped  that  the 
material  may  be  used  as  a  practical  reference  for  those  gaining  their 
first  exposure  to  the  subject,  although  key  references  are  given  for 
further  research  into  specific  points.  Any  errors  of  omission  most 
likely  reflect  the  author's  limited  exposure  to  the  field  through 
his  application  of  the  method  to  a  few  particular  research  problems. 
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ABSTRACT 


A  short  presentation  is  made  on  the  basic  elements  of  Power 
Spectral  Analysis  with  emphasis  on  the  Blackman-Tukey  method.  Short 
discussions  are  included  on  the  topics  of  pre-whitening,  frequency 
and  spectral  windows,  and  statistical  reliability.  Examples  are 
Included  whenever  possible,  and  a  Fortran  subroutine  for  calculating 
a  power  spectrum  is  presented. 
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I.  BRIEF  BACKGROUND 

A.  Historical  Roots 

Power  spectral  analysis  is  one  small  area  of  the  much  broader 
field  of  Communications  Theory.  This  broader  field  is  an  indispens¬ 
able  part  of  communications  engineering  and  provides  the  theoretical 
foundations  for  the  design  and  analysis  of  much  of  the  advanced  engi¬ 
neering  found  in  modern  communications  systems  (Lee,  i960).  The 
theory  is  basically  a  statistical  theory  in  which  the  central  idea  Is 
that  noise  and  messages  are  considered  to  be  random  phenomena.  Prob¬ 
ability  theory  is  therefore  incorporated  into  the  very  foundation  of 
the  theory  and  is  an  Integral  part  of  it. 

The  basis  of  the  theory  finds  its  roots  in  statistical  mechanics. 
The  equivalence  of  time  and  ensemble  averages,  first  assumed  by  Gibbs 
and  later  stated  more  precisely  by  Maxwell  in  his  ergodlc  hypothesis , 
is  the  starting  point  for  statistical  communications  theory,  from 
this  and  the  quaai-ergodlc  hypothesis  are  derived  the  formal  proofs 
necessary  for  the  logical  development  of  the  subject. 

B,  Ergodiclty  and  Stationarity 

Two  conditions  necessary  for  the  development  of  the  theory  are 
imposed  on  the  random  ensembles  of  data  which  we  wish  to  power  spec¬ 
trum  analyze.  We  shall  merely  state  them  as  being  the  foundation  for 
the  development  of  the  theory,  with  proofs  and  implications  to  be 


found  elsewhere. 
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(l)  Ergodicity.  The  ergodic  theorem  may  be  stated  as  "in  a 
stationary  ensemble  of  random  functions  having  a  continuous  range  of 
possible  values  the  amplitudes  of  an  ensemble  member  will  come  infin¬ 
itely  close  to  every  point  of  the  continuous  range  of  possible  values 
if  given  an  infinite  amount  of  time."  This  theorem  allows  the  replace 
ment  of  ensemble  averages  with  time  averages,  and  is  the  basis  for  the 
formal  analysis  of  communications  theory. 

(?)  Stationarity.  If  the  amplitude  probability  density  of  an 
ensemble  is  time  independent,  the  ensemble  is  said  to  be  stationary. 

In  practical  terms  related  to  power  spectral  analysis  this  means  that 
the  power  spectrum  of  a  finite  data  set  is  time  independent. 

In  addition  to  the  two  above  conditions,  it  is  assumed  that 
the  random  process  is  Gaussian  or  nearly  Gaussian  in  character,  that 
is,  the  probability  distribution  of  the  elements  in  an  ensemble  is 
Gaussian  or  nearly  so.  Blackman  and  Tukey  (1958)  show  that  for  an 
infinite  data  set  a  Gaussian  assumption  yields  exact  results,  and 
rather  good  approximations  otherwise. 

C,  Notation 

With  these  preliminaries  stated,  we  now  give  the  notation  to 
be  used  throughout  the  rest  of  this  section.  Fourier  transforms, 
correlations  (auto-  and  cross-),  and  convolutions  are  used  rather 
frequently  in  power  spectrum  work.  Therefore  to  reduce  the  complexity 
of  the  equations  in  the  following  disoussion  a  simplified  notation 
will  be  adopted  as  follows: 


1.  Fourier  Transform 


Let  f(t)  be  specified  on  the  interval  (-«,  ®).  ltoen  define 
the  Fourier  transform 

F(o»)  *  —  j  *  f(t)  e’iut  dt  . 

m  - 

We  shall  also  use  the  Fourier  transform  operator  e.g., 

*[f(t)]  =  F («) 


and  Inverse  Fourier  transform  operator 


3rlCf(t)i  -  -i-  r*  f(t)  ^  dt 

/2n  * 


2.  Correlation  Functions 

Let  f^(t)  and  f?(t)  be  specified  on  the  interval  (-•,  •).  We 
then  use  the  following  notation: 
a.  Convolution.  Define 


♦12(t)  2  11m  fL(t)  tp(r  -  t)  dt 

=  fL(t)  *  f2(t) 


where  *  denotes  convolution 


b .  Cross  Correlation.  Define 


v(t)  '  fi(t)  +  t)  dt 

3  fL(t)  *  f2(-t)  . 

c.  Autocorrelation.  Cross  correlation  of  a  function  with 
itself.  Define 

Vll^  =  X-t/p  V*  +  T)  dt 

-  fx(t)  *  fx(-t) 

=  fx(t)  *  fx(t)  , 

the  last  step  being  a  result  of  the  symmetry  of  the  autocorrelation 
operation. 

In  general,  then,  small  letters  will  denote  functions  of  time 
and  capital  letters  Fourier  transforms  (functions  of  frequency)  of 
the  corresponding  functions  of  tire.  Letters  written  in  script  will 
be  used  to  denote  operators.  Additionally,  MLP  will  be  used  as  an 
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abbreviation  for  mean-lagged-product,  and  FFT  for  fast-Fourier- 
transform.  The  term  "power-spectral-density"  (PSD)  will  also  be 
used  synonymously  with  "power  spectrum". 


6 


II.  THE  POWER  SPECTRUM  FOR  THE  CONTINUOUS  CASE 


A  Fourier  series  is  one  set  of  orthogonal  functions  that  may 
be  used  to  expand  a  well-behaved  function  on  the  interval  (-•,  •). 
If  we  consider  a  time  series  f(t),  we  can  represent  it  by 

f(t)  =  ^1CFC«))] 


where 


F(u>)  =  7[f(t)]  . 

The  power  spectrum  of  a  function  is  defined  as  the  absolute  value 
squared  of  the  Fourier  transform  of  the  function.  If  P(u>)  is  the 
power  spectrum,  then 

PM  «  l>ff(t)]|2  »  |fWIs  . 

But  because  of  the  convolution  theorem,  which  states  that  for  two 
functions  f^  and  fp, 

?Ff  *  f  1  =  F  •  F 
1  1  2J  1  e2 
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By  setting  =  f?  we  have 

,ffl  *  fl]  =  Fi  *  Fi  * 

Now  f^  *  is  the  autocorrelation  function  for  f^. 
tion  function  is  always  an  even  function,  i.e., 

fx(t)  *  fx(-t)  *  fx(t)  *  fj_(t) 

so  that  there  are  no  sine  components  in  its  Fourier 
transform  is  therefore  real  and 


»[fx  *  fxl 


»fvu]  - 


T11 


so  that 


P(«)  -  |F(u>)|2 


(2) 

The  autocorre  la- 


transform.  The 


(5) 


00 
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Therr  '’ore 


P(tu)  ■  $(u>)  = 


(5) 


or 


» 

*|f(u>)|?  =  ^(ou)  .  (6) 

This  equality  is  known  as  the  Wiener  Theorem  and  states  that 
the  power  spectrum  of  a  time  series  is  equal  to 

(a)  the  Fourier  transform  squared  of  the  function,  or 

(b)  the  Fourier  transform  of  the  autocorrelation  function. 

It  may  be  noted  that  in  (a)  the  Fourier  transform  utilizes  both 
sine  and  cosine  terms,  while  in  (b)  the  Fourier  transform  of  the 
autocorrelation  function  has  only  non-zero  cosine  terms,  that  is, 

=  *r<pu(tn 


j.. 


^li (t)  cos  ^ 


P  •  ^ 

-  J  <c,,(t)  co  ut  dt 

0  11 


(7) 
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The  first  method  of  computing  the  power  spectrum  is  the  most 
direct  and  straightforward  way.  The  fast- Fourier- transform  (FFT) 
calculation  utilizes  the  definition  of  the  power  spectrum  directly 
to  compute  the  power  spectrum.  The  second  method  is  known  as  the 
mean-lagged-product  (MLP)  method,  and  is  the  one  utilized  by  Blackman 
and  Tukey  (1958)  to  calculate  power  spectra. 
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III.  FINITE  DATA  SETS 


A.  Effect  of  Data  Truncation  on  the  Power  Spectrum 
The  definition  of  the  power  spectrum  was  made  for  infinitely 
long,  continuous  data  sets,  fractical  data  analysis,  however,  re¬ 
quires  the  use  of  considerably  less  data.  To  determine  the  effects 
on  the  power  spectrum  resulting  from  the  truncation  of  a  data  set 
to  finite  lengths,  consider  the  function  f(t)  defined  on  the  inter¬ 
val  (-•,  •)  as  shown  in  Figure  1. 


f(t) 


- . - ♦ 

Figure  1 


If  we  truncate  f(t)  by  multiplying  by  a  data  window  g(t)  such  that 


g(t) 


S 


It  |  s  T/2 

it  I  >  T/2 


the  truncated  data  set  becomes 


h(t)  -  f(t)  g(t) 
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as  shown  in  Figure  2. 


The  function  h(t)  then  represents  the  truncated  data  set  available 
from  which  the  power  spectrum  is  to  be  calculated.  The  power  spectrum 
P  (u>)  of  h(t),  representing  the  apparent  power  spectrum  of  f(t),  is 
then 

P  (<*>)  -  l*(f  *  g)l? 

.  |f  *  g12 

-  If |2  *  |ol2  ,  (8) 

the  last  step  resulting  from  the  fact  that  g  is  presumed  to  be  an 
even  function.  Thus  the  true  power  spectrum 

W”’  -  |r|2 
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is  modified  by  convolution  with  |G|'  ,  the  Fourier  transform  squared 
of  the  data  window.  G  is  called  the  frequency  window  by  Blackman  and 
Tukey  (1958).  For  the  data  set  shown  in  Figure  A-?,  G  is  the  sine 
function,  shown  in  Figure  5. 


Convolution  by  the  frequency  windtv  causes  a  certain  degree  of 
smoothing  in  the  calculated  power  spectrum  and  a  small  amount  of 
leakage  via  the  side  lobes  from  nearby  frequency  bands  into  the  fre¬ 
quency  band  of  interest.  If  the  confuted  power  spectrum  is  relatively 
flat,  i.e.,  has  a  dynamic  range  of  less  than  2  or  5  orders  of  magni¬ 
tude,  this  smearing  or  leakage  causes  little  or  no  problem,  for  the 
amplitude  of  the  largest  (first)  side  lobe  of  the  sine  function  is 
only  on  the  order  of  several  percent  of  that  of  the  main  lobe.  If, 
however,  there  are  large,  well-defined  peaks  in  the  power  spectrum 
such  large  peaks  act  as  a  first  approximation  delta  function.  When 
convolved  with  the  frequency  window  they  act  to  reproduce  the  win¬ 
dow,  producing  spurious  peaks  in  the  power  spectrum  corresponding  to 
the  side  lobes  of  the  frequency  window.  These  spurious  peaks  may 
be  mistakenly  identified  as  structural  details  of  the  true  power 
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spectrum  when  in  reality  they  are  merely  artifacts  created  by  trunca¬ 
tion  of  the  original  data  set. 

B,  Reducing  Frequency  Window  Leakage 
Because  of  the  frequency  window  side  lobe  leakage  associated 
with  data  function  truncation,  one  is  properly  concerned  with  means 
of  reducing  or  minimizing  the  undesirable  effects  of  such  truncation. 
Several  methods  are  available  for  doing  this,  each  depending  somewhat 
on  the  particulars  of  how  the  power  spectrum  is  computed.  Three 
commonly  used  methods  will  be  described  in  the  following  sections, 
each  being  intended  to  illustrate  the  basic  features  of  and  rationale 
behind  each  procedure. 


1.  Data  Tapering 

One  of  the  obvious  methods  of  reducing  the  side  lobe  leakage 
is  to  choose  the  data  window  g(t)  such  that  its  Fourier  transform 
G(uj)  has  either  no  aide  lobes  at  all  or  side  lobes  that  are  small 
compared  to  the  main  lobe.  Examples  of  the  former  are: 


g(t) 


.(|t|/Sto)2 

e 


(Gaussian) 


where  t  is  some  scale  factor.  This  function  has  a  Fourier  transform 
o 

equal  to  a  Gaussian;  and 


sin  2tt  Af  t 


g(t) 


Pn  Af  t 


(sine  function) 


which  has  a  Fourier  transform  equal  to  the  box-car  function. 

However,  for  the  side  lobes  to  be  eliminated  altogether  it  is 
necessary  to  extend  these  data  windows  to  ±  «*.  If  they  are  truncated 
at  some  finite  value  (as  they  must  be  in  practice),  side  lobes  are 
again  introduced,  though  they  will  be  smaller  than  the  case  where 
g(t)  is  the  box-car  function. 

A  simpler  method  is  to  taper  the  data  set  at  the  ends  using 
data  windows  indicated  in  the  following  sketches: 


or 


Although  side  lobes  are  still  present  in  the  frequency  windows 
for  each  of  the  above  cases,  they  will  be  smaller  than  for  the  case 
where  the  simple  box-car  data  window  was  used.  The  smaller  side  lobes 
are  a  result  of  replacing  the  abrupt  discontinuities  of  the  original 
box-car  data  window  with  more  gently  sloping  functions. 
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Data  tapering  is  most  often  performed  vhen  the  FFT  method  is 
used  to  compute  the  power  spectrum.  For  a  more  conplete  discussion 
of  this  topic  see  Ehochson  and  Otnes  (1968). 

?.  Tapering  the  Autocorrelation  Function 

If  the  MLP  method  of  computing  the  PSD  is  used  the  problem  of 
side  lobe  leakage  is  not  as  simple  as  it  was  in  the  case  of  computing 
the  PSD  directly  from  its  definition.  An  additional  factor  for  con¬ 
sideration  enters  when  not  only  the  original  data  set  must  be  trun¬ 
cated,  but  so  also  the  autocorrelation  function.  When  an  MLP  cal¬ 
culation  is  made  it  becomes  advisable,  for  reasons  to  be  discussed 
in  connection  with  statistical  reliability,  to  truncate  the  auto¬ 
correlation  function  at  a  maximum  lag  of  not  more  than  10—  20f,  of 
the  length  of  the  data  set  (Blackman  and  Tukey,  1958).  In  order  to 
see  what  effect  this  second  truncation  has,  consider  the  original 
data  set  f(t)  properly  truncated  with  a  data  window  g(t).  The  auto¬ 
correlation  function  then  becomes 

*u(t)  *  [f(t)  •  g(t)]  *  [f(t)  •  gft)] 

or 

’ll  "  *  g)  *  (f  ’  g) 


4 

1 
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Now  cp^  represents  the  entire  autocorrelation  function  avail- 
able  after  f  is  truncated  by  g.  To  truncate  at  a  lag  of  10— 20fj6 
of  the  length  of  the  data  set  therefore  involves  specifying  another 
function  q(r),  called  the  lag  window ,  by  which  is  multiplied 
in  order  to  effect  truncation.  Let  rp^  represent  the  truncated 
autocorrelation  function.  Then 

®11  “  *11  *  q  "  Kf  *  g)  *  (f  *  8)1  *  <3 

Hie  apparent  power  spectrum  is  the  Fourier  transform  of  the 
truncated  autocorrelation  function  so  that 

Pap(u>)  =  *{[(f  •  g)  *  (f  •  g)]  •  q} 

-  >[(f  •  g)  *  (f  •  g)l  *  *(a) 

*  Wf  •  g)  *  *(f  •  g)l  *  y(q) 

•  r(F  *  G)  •  (F  *  0)1  *  Q  ,  (9) 

where  Q  (the  so-called  spectral  window)  is  the  Fourier  transform  of 
q,  the  lag  window.  Now  g  and  q  are  normally  chosen  to  be  even 
functions,  so  that 
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PU)  =  ( |FlP  *  jc  1P)  *  Q 

ap 

=  fptrue(«>)  *  iG((u)  |?]  *  Q  .  (10) 

We  see  that,  as  In  section  III-A,  the  true  power  spectrum  has  been 
altered  by  convolution  with  |g(u>)  |?,  and  additionally  by  convolution 
with  Q.  If  q(T)  was  the  box-car  function,  then  Q(t»)  is  the  familiar 
sine  function,  not  squared.  It  is  this  last  convolution  that  can 
lead  to  negative  values  in  the  apparent  power  spectrum  even  though 
negative  values  are  theoretically  impossible  in  a  power  spectrum. 
They  are,  in  this  case,  merely  an  artifact  due  to  truncation. 

To  remove  the  possibility  of  negative  values  for  the  power 
spectral  estimates,  as  well  as  making  the  side  lobes  of  the  window 
Q  smaller,  it  is  once  again  advisable  to  tailor  the  shape  of 
the  lag  window  o(t).  It  is  to  be  tailored  in  such  a  way  that  the 
side  lobes  that  remain  are  small  in  comparison  to  the  main  lobe,  end 
damp  out  quickly  with  increasing  distance  from  the  main  lobe.  Two 
lag  windows  that  have  been  found  to  do  this  effectively  are: 


(1)  c^t)  .  ill  ♦  CO,  f 

\  m 


•  T  i  t  s  T 


m 


m 


where  T  is  the  greatest  lag  used  in  the  autocorrelation  function. 

IQ 
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The  use  of  this  window  is  called  "hanning".  A  second,  more  conmonly 
used  lag  window  is 

(2)  q?(r)  »  0.54  +  0.46  cos  ^  ,  -  Tm  <  T  5  Tm  * 

m 

Use  of  this  window  is  called  "hamming". 

The  two  functions  q^  and  q^,  along  with  their  associated 
transforms  Q1  and  are  shown  in  Figure  4. 


Figure  4 

(After  Blackman  and  Tukey,  1959.) 
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Although  in  principle  q(T)  would  normally  be  applied  to  the 
autocorrelation  function  prior  to  transforming,  the  interchange¬ 
ability  of  the  integrations  associated  with  the  transform  and  the 
convolution  makes  it  possible  to  convolve  with  Q((l))  after  the 
transform  has  been  completed.  This  has  particular  advantages  when 
the  data  are  in  digital  form.  In  the  case  of  hamming  or  harming, 
the  convolution  will  take  the  form  of  a  5-point  smoothing  formula 
easily  applied  to  the  power  spectrum  estimates  calculated  by  trans¬ 
forming  the  autocorrelation  function.  More  will  be  said  about  this 
form  of  smoothing  in  section  V  when  a  specific  power  spectrum  exanple 
is  given. 

It  might  be  added  that  when  the  MLP  method  for  computing  PSD's 

is  used  it  is  usually  unnecessary  to  be  overly  concerned  about  the 

2 

effects  of  the  convolution  of  |Q(f„)l  with  Ptrue(«»)  (see  Eq.  (10)]. 

2 

Since  G((U)  is  squared,  whereas  Q((«)  is  not,  the  side  lobes  of  |G(,„)| 
will  be  unimportant  compared  to  those  of  Q(u>).  The  half-width  of  the 
major  lobe  of  G((u)  is  also  small  (~  1  order  of  magnitude  smaller) 
compared  to  that  of  Q(«,)  since  the  original  data  set  is  ~  10  times 
longer  than  q (t)  Therefore  the  effects  of  Q  will  far  outweigh  those 
of  G,  and  it  is  those  effects  that  hamming  or  banning  are  designed 
to  offset.  It  is  accordingly  not  necessary  to  taper  the  original  data 
set  (or  pre -whiten;  see  next  section)  except  in  cases  of  extremely 
discontinuous  spectra,  or  when  extreme  care  need  be  taken  to  assure 
the  validity  of  the  results. 
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3.  Pre -whitening 

As  was  mentioned  previously,  if  there  are  large,  well-defined 
peaks  in  the  power  spectrum  such  peaks  can  produce  spurious  detail 
in  the  power  spectrum  due  to  frequency  or  spectral  window  side  lobe 
leakage.  It  is  for  this  reason  that  methods  were  sought  to  reduce 
the  side  lobes.  An  alternate  method  that  can  be  used  on  occasion  is 
known  as  "pre -whitening".  This  procedure  involves  flattening  the 
power  spectrum  prior  to  calculation  by  passing  the  data  through  a 
filter  with  a  known  power  transfer  function  in  order  to  eliminate 

large  peaks  and  discontinuities.  With  the  peaks  and  gross  discontin- 

2  2 

uities  thus  removed,  the  effective  convolution  of  |F|  with  |G|  to 
yield  Pap(m)  will  not  act  to  reproduce  the  side  lobes  of  |G|  as 
would  have  been  the  case  had  the  peaks  not  been  removed.  Once  the 
power  spectrum  is  computed  the  inverse  of  the  pre -whitening,  the 
so-called  "post-darkening",  is  applied  to  the  PSD  estimates  to 
complete  the  calculation. 

In  order  for  pre-whitening  to  be  used  effectively  some  prior 
knowledge  of  the  expected  shape  of  the  spectrum  to  be  flattened  mu$t 
be  available.  This  is  necessary  in  order  to  design  a  filter  with 
the  proper  power  transfer  function.  An  example  of  how  pre-whitening 
may  be  used  is  the  following:  Suppose  wo  have  a  data  set  for  which 
we  know  approximately  the  shape  of  the  power  spectrum,  and  suppose  it 
has  a  very  large  low  (zero)  frequency  component,  as  in  Figure  5. 


Figure  5 


If  the  power  spectrum  is  computed  directly  from  the  data  with¬ 
out  some  measures  being  taken  to  compensate  for  possible  side  lobe 

leakage,  the  apparent  power  spectrum  P  (ju)  would  look  something  like 

ap 

Figure  6. 


Figure  6 


p? 


The  small  scale  structure  may  well  be  that  solely  due  to  side 
lobe  leakage,  though  in  the  case  where  confidence  in  the  results  is 
low  this  structure  can  be  masked  by  statistical  noise. 

Suppose  now  that  the  original  data  f ft)  is  pre-whitened  by 
convolution  with  a  smoothing  function  sft).  Then  the  data  set  h(t) 
available  from  which  the  PSD  calculation  is  made  will  be 

hft)  =  f(t)  *  sft) 

or  more  simply, 

h  =  f  *  s 


Then 


P  («)  =  l*fh)i2  =  i?(f  *  s)|2 
ap 


or 


ap 


=  F 


Ptrue((u) 


(11) 


since  s(t)  is  an  even  function. 

Thus,  if  the  original  data  set  is  convolved  with  a  smoothing 
function  s(t),  the  true  PSD  is  modified  by  the  direct  product  of  the 
square  of  the  transform  of  s(t).  |s|2  is  called  the  power  transfer 

function.  To  compensate  for  this  at  the  end  of  the  calculation  we  mul¬ 
tiply  by  the  inverse  of  |s(u>)  |?  to  restore  the  proper  shape  to  P(u>). 

In  the  above  example  suppose  we  choose  a  smoothing  function 

It  -  tQ|  < 

It  -  t  |  *  ?At 
o 

where  At  is  some  scale  factor.  This  smoothing  function  may  be  passed 
over  the  data  as  many  times  as  one  desires,  more  smoothing  being 
accomplished  with  each  pass  and  more  of  the  high  frequency  components 
being  suppressed. 

If  the  data  is  in  digital  form  the  smoothing  function  above 
takes  the  form  such  that  if  f*(t)  is  the  smoothed  data, 

fi(°  *  ?  fi(t)  +  l  rfi-lft)  +  fi+i(t)1  *  (1?) 

Holloway  (1958)  showed  that  for  n  successive  passes  of  the  above 
elementary  filter  function  through  the  data,  the  power  transfer  func¬ 
tion  |s|?  becomes 


1  +  cos 


n(t  -  tQ) 


PAt 


i  (t ) 


|S |^  =  cos  ^(nf  At) 


2k 


where  At  is  the  data  sample  spacing. 

The  above  filter  function  may  be  used  to  effectively  remove 
the  low-frequency  component  in  the  present  example.  If  the  smoothed 
data  (low-pass  filtered)  is  subtracted  from  the  original  unsmoothed 
data,  the  difference  will  represent  the  original  data  filtered  by  a 
high-pass  filter  with  a  power  transfer  function  equal  to  the  comple¬ 
ment  of  the  original  transfer  function.  Suppose  the  original  data 
is  filtered  using  this  method.  The  smoothing  function  applied  to  the 
original  data  set  f(t)  will  have  a  power  transfer  function  similar  to 
Figure  7. 


When  this  smoothed  data  is  subtracted  from  the  original  f(t), 
the  resultant  power  transfer  function  will  be  the  complement  of  the 


function  shown  above. 
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Applied  to  our  example  of  data  with  a  large  low-frequency 
component,  this  will  effectively  flatten  the  curve  to  minimize  the 
possibility  of  side  lobe  leakage.  After  the  PSD  has  been  computed 
post-darkening  is  achieved  by  multiplying  by  the  inverse  of  the  func¬ 
tion  shown  in  Figure  8  to  complete  the  calculation. 

For  a  more  complete  discussion  of  the  topic  of  pre-whitening 
and  digital  filtering  techniques  see  BlaeOTan  and  Tukey  (1958)  and 
Qiochson  and  Otnes  (19^8). 
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IV.  DIGITIZED  DATA 

The  discussion  of  power  spectra  has  been  general  up  to  this 
point,  with  only  an  occasional  reference  to  specifies.  Since  we  are 
here  primarily  interested  in  data  that  appears  in  discrete,  digital 
form,  it  is  appropriate  to  specialize  to  that  case.  We  shall  from 
this  point  forward  consider  only  the  power  spectra  of  data  consisting 
of  discrete,  eaui-spaced  samples.  The  discrete  forms  of  the  general 
equations  of  section  II  will  be  given  for  both  the  MLP  and  FFT  cal¬ 
culations.  The  statistical  reliability  of  PSD  estimates  will  be 
discussed  briefly  for  each  of  the  two  metiods,  and  several  considera¬ 
tions  helpful  for  planning  will  be  mentioned.  The  averaging  of 
several  power  spectra  is  also  mentioned. 

A.  Discrete  Forms  of  Pelevart  Equations 

1.  MLP  Method 

The  steps  required  for  calculation  of  a  power  spectrum  using 
the  MLP  method  may  be  summarized  from  the  above  discussion  to  be  the 
following : 

(l)  Pre-whitening.  If  the  power  rpectrum  is  known  to  contain 
large  peaks  or  discontinuities,  the  raw  data  should  be  pre-whltened 
by  use  of  appropriate  digital  filters  (high-pass,  low-pass,  band-pass, 
or  combinations  of  these).  The  necessity  for  this  procedure  is 
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somewhat  open  to  interpretation  according  to  one's  perception  of 
what  constitutes  a  large  peak  or  discontinuity.  My  own  limited 
experience  over  several  years  is  that  for  a  power  spectrum  with  a 
dynamic  range  of  >  2—3  orders  of  magnitude,  pre-whitening  to 
flatten  the  spectrum  prior  to  its  computation  is  advisable.  The 
example  given  in  section  III  of  a  method  of  constructing  a  high- 
pass  filter  was  used  successfully  by  the  author  (M.S.  thesis,  1973). 
Since,  however,  the  actual  method  used  to  pre-whiten  will  depend 
on  the  details  of  the  various  power  spectra  encountered  in  practice, 
the  reader  is  referred  to  chapter  3  of  Enochson  and  Otnes  (1968) 
for  a  thorough  discussion  of  recursive,  non-recursive,  and  second- 
and  higher-order  filters  and  their  application  to  time  series. 

(?)  Normalizing  the  Data.  Although  not  mentioned  ii  the 
previous  sections,  it  is  advisable  in  practice  to  normalize  the  data 
to  zero  mean  and  unit  standard  deviation  before  calculation  of  the 
PSD.  Since  the  calculation  of  the  mean-lagged-products  (autocorrela¬ 
tion  function)  involves  the  sum  of  many  products,  it  is  easy  to  termi¬ 
nate  a  computer  calculation  of  a  PSD  prematurely  due  to  overflow. 
Subtracting  the  mean  from  the  data  removes  only  the  zero-frequency 
cosine  term  from  the  PSD,  and  can  be  included  in  the  calculation  after 
the  remaining  following  steps  are  completed.  Dividing  the  (pre¬ 
whitened)  zero-mean  data  set  by  o  reduces  the  data  to  unit  standard 
deviation.  Correct  absolute  units  can  be  restored  to  the  final  com- 

r 

puted  PSD  if  desired  by  multiplying  by  d  . 
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(7>)  Calculation  of  the  autocorrelation  function.  If  tp. 

^  j 

denotes  the  j'th  value  of  the  autocorrelation  function,  the  discrete 
form  for  the  autocorrelation  function  of  the  data  set  f(t)  with  data 
sample  spacing  At  containing  n  discrete  values  is 

*3  =  rH  Jjj  fi  '  fi+j  ’  j  =  0,  1,  2,  ...  m  ,  (13) 

where  m  is  the  maximum  number  of  lags  in  the  autocorrelation  function. 
As  mentioned  in  section  II,  m  will  generally  be  limited  to  10— 20J& 
that  of  n. 

(4)  Fourier  transform  of  the  autocorrelation  function.  Since 
the  autocorrelation  function  is  an  even  function,  the  discrete  form 
for  the  Fourier  transform  becomes  the  discrete  finite  cosine  series. 
Applying  this  to  the  sequence  v  .  fa. ,  . . . ,  v  we  obtain  as  raw 
estimates  for  the  PSD 


At 


(14) 


i  =  0,  1,  ...  m 


where  At  is  the  sample  spacing. 
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(5)  Smoothing  the  raw  estimates.  The  raw  estimates  calculated 
in  step  (4)  correspond  to  Eq.  (10) 

%<”'>  '  rPtrue(»)  *  *  «  - 

with  a  box-car  lag  window  leading  to  a  sine  function  spectral  window. 
What  remains  is  convolution  with  a  suitably  tailored  spectral  window 
Q  with  small  side  lobes  to  complete  the  calculation.  A  commonly 
used  window  is  the  hamming  window,  the  digital  form  of  which  becomes 

?i  =  0.54  P1  +  0.2^5  (Pi+1  +  PU1)  .  (15) 

(C)  Post -darkening.  If  the  raw  data  were  pre-whitened  in 
step  (l),  the  last  step  in  the  calculation  of  the  PSD  will  be  to 
restore  the  true  shape  of  the  power  spectrum  by  post-darkening.  This 
process  involves  multiplying  the  smoothed  estimates  of  step  (5)  by 
the  inverse  of  the  pre-whitening  filter  power  transfer  function. 

The  frequency  resolution  Af  of  the  resultant  spectrum  will  be 


Af 


1 

2m  At 


where 

m  =  the  maximum  number  of  lags  in  the  autocorrelation  function 
end 


At  =  data  set,  sample  period. 
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The  maximum  frequency  of  the  computed  PSD  will  be 


f 

max 


=  irAf  = 


1 

2At 


(16) 


in  accordance  with  the  sampling  theorem. 

?.  FFT  Method 

The  fast -Fourier- transform  technique  of  computing  the  PSD  of 
a  time  series  was  a  major  improvement  over  the  MLP  method  in  terms  of 
the  actual  computer  time  taken  for  the  calculation.  If  N  is  the  num¬ 
ber  of  points  in  a  time  series,  the  total  computer  time  needed  for 

p 

an  MLP  calculation  is  roughly  proportional  to  Nr ,  where  for  the  FFT' 
method  it  goes  approximately  as  N(log  N).  The  advantages  of  utilizing 
the  FFT  method  whenever  many  points  are  being  power  spectrum  analyzed 
therefore  lie  on  the  side  of  efficiency  rather  than  any  fundamental 
superiority  of  the  method  over  that  of  the  MLP.  As  a  technique  of 
computing  PSD's,  it  is  quickly  supplanting  the  MLP  method  and  is 
therefore  worth  studying.  An  extensive  discussion  of  the  FFT  com¬ 
puted  one-dimensional  PSD  by  Brault  and  White  (1971 )  provides  a 
thorough  introduction  to  FFT  methods  and  their  application  to  astro¬ 
nomical  problems. 

Although  no  attempt  will  be  made  here  to  outline  in  detail  the 
steps  required  for  FFT  calculation  of  a  PSD,  the  basic  steps  are  as 


follows : 
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(l)  Pre-whitening.  The  necessity  for  pre-whitening  the  data 
is  less  strong  when  the  FFT  method  is  used  than  it  is  for  the  MLP 
method.  The  reason  for  this  lies  in  the  fact  that  the  FFT  calculated 
PSD  takej  the  form 


P 

ap 


(») 


p,  * 
true 


|g(«»)  l? 


9 


whereas  the  MLP  method  yielded 

P  (w)  -  rPtrue  *  |G(a>)|2]  *  Q(tti) 


The  FFT  form  does  not  involve  a  convolution  with  Q,  so  that 
one  need  only  be  concerned  with  the  shape  of  G(u>).  If  g(t)  is 
properly  specified  the  side  lobes  of  |G(u,)  |  can  be  kept  small  enough 
so  that  very  little  leakage  is  present  [s?e  step  (3)  below].  Only 
in  the  case  where  there  are  exceptionally  large  peaks  or  discontin¬ 
uities  in  the  power  spectrum  should  pre -whitening  be  necessary.  In 
general,  it  may  be  said  that  in  most  cases  this  step  will  not  be 
necessary  at  all. 

(2)  Normalizing  the  data.  It  is  advisable  to  normalize  the 
data  to  zero-mean  and  unit  standard  deviation,  for  the  same  reasons 
as  given  for  the  MLP  method  above. 

(3)  Data  window  correction.  As  discussed  in  section  III-A, 
truncation  of  a  data  set  can  produce  spurious  detail  in  the  computed 
PSD.  It  is  therefore  necessary  to  choose  a  data  window  g(t)  by  which 
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to  multiply  the  time  series  such  that  the  frequency  window  G(iu)  will 
have  small  side  lobes.  A  window  currently  in  common  use  is  the  ten- 
percent  cosine  bell,  the  application  of  which  is  called  coswlnding. 

For  the  generalized  formula  for  this  data  window  see  Brault  and  White 
(1971,  Equation  13). 

(4)  Appending  zeros  to  the  data.  The  most  efficient  versions 
of  the  FFT  require  the  data  to  consist  of  a  number  of  points  equal 

to  a  power  of  two,  although  some  generalized  versions  will  work  with 
an  arbitrary  number  of  points.  If  a  version  is  used  requiring  some 
specific  number  of  points,  the  original  normalized  data  set  to  be 
transformed  must  have  zeros  appended  to  it  to  bring  the  total  points 
up  to  the  specified  number. 

(5)  Tr^£fqrrinB.the.dat?.  The  normalised  date  aet  f,  !e 
transformed  using  the  forward  discrete  Fourier  series  transform  to 
yield  Fourier  coefficients 


1.  +  ib.  =  —  Y  f.  exp  i  -jk  — — =- 
J  j  n^  .Cj  k  ^  0  n  -  1 

O  K*0  ,  O 


V1 


,  J  =  0,  1,  • .  • ,  nQ  -  1 


(17) 


where 


At  =  data  sample  spacing, 

nQ  =  number  of  data  points  in  the  normalized  data  set,  and 
i  =  /Tl  . 

(6)  Computing  the  raw  spectrum.  The  raw  PSD  is  Just  the  ab¬ 


solute  value  squared  of  the  Fourier  transform,  adjusted  to  compensate 
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for  the  appendage  of  zeros  to  the  original  data  set. 


2 

no  /  2 
n  ^ai 


+  bj) 


(18) 


where  n  *  number  of  points  in  the  original  data  set. 

(7)  Smoothing  the  raw  estimates.  The  raw  estimates  obtained 
in  step  (6),  although  theoretically  correct,  may  be  statistically 
unreliable.  It  is  possible  to  increase  the  statistical  reliability  by 
smoothing  (convolving)  the  estimates  with  a  suitable  smoothing  func¬ 
tion.  For  an  excellent  discussion  on  the  subject  of  smoothing  raw  FFT 
estimates  see  Edmonds  and  Webb  (197?). 

The  frequency  resolution  (PSD  estimate  spacing)  of  the  FFT 
calculated  PSD  will  be 


Af 


1 

n  At 
o 


(19) 


For  specifics  on  the  programming  of  the  FFT  algorithm,  including 
Fortran  indexing  peculiarities,  one  should  consult  'Special  Issue 
on  FFT  and  Its  Applications  to  Digital  Filtering  and  Spectral  Analysis', 
IEEE  Trans.  AJM£,  No.  2  (1967).  Another  special  issue  on  the  same 
topic,  IEEE  Trans.  AU^IJ,  No.  ?  (19^9) »  gives  an  extensive  biblio¬ 
graphy. 


B.  The  Statistical  Reliability  of  the  PSD  Estimates 


1.  Confidence  Limits 

The  calculation  of  the  PSD  from  a  finite  number  of  discrete 
data  points  can  be  expected  to  have  associated  with  it  a  set  of  confi¬ 
dence  limits  reflecting  the  presumed  statistical  nature  of  the  ori¬ 
ginal  time  series.  In  theory,  if  a  PSD  calculation  could  be  made  from 
an  infinite  number  of  data  points,  one  would  have  absolute  confidence 
in  the  results.  In  practice  one  must  settle  for  a  finite  number  of 
points  ftom  which  to  make  the  calculation.  Each  point  in  the  resul¬ 
tant  power  spectrum  will  have  a  set  of  "confidence  limits"  assigned 
to  it  reflecting  the  confidence  (in  a  statistical  sense)  that  the  re¬ 
sulting  power  spectrum  is  not  due  solely  to  randomly  distributed 
data  points. 

The  statistical  accuracy  of  the  computed  PSD  is  estimated  in 
terms  of  "equivalent  degrees  of  freedom",  from  which  the  confidence 
limits  are  calculated.  The  degrees  of  freedom  may  be  thought  of  as 
representing  the  number  of  estimates  of  power  in  the  frequency  inter¬ 
val  Af,  the  frequency  resolution  of  the  computed  PSD.  A  measure  of 
the  probability  that  an  estimate  falls  within  an  upper  and  lower 
bound,  the  ratio  of  which  is  designated  the  confidence  factor  fc>  is 
given  by 


f 


c 


10b/aa/ro 


exp 


2.3  b 
,10/k  -  1 


(?0) 


where 


f  =  confidence  factor  for  a  given  confidence, 
k  =  degrees  of  freedom,  and 

b  =  factor  depending  on  the  desired  confidence,  given  as 
follows : 


confidence 

50$  80$  90$  96$  98$ 

b  [~~8  16  20  25  29 

This  is  an  approximation,  but  for  k  ^  4  is  very  close  to  more  exact 

calculations  based  upon  chi-square  tables  (Edmonds,  1966). 

The  confidence  factor  ffi  may  be  used  as  error  bars  for  the 

computed  PSD,  being  positioned  on  each  data  point  in  such  a  way 

that  the  point  flails  on  the  geometric  mean  of  the  upper  and  lower 

limits  of  f  .  Another  way  of  stating  this  is  that  the  upper  limit 

for  the  confidence  limits  is  equal  to  the  PSD  estimate  times  the 

square  root  of  f  ;  the  lower  limit  equals  the  estimate  divided  by  the 

square  root  of  f  .  If  the  PSD  is  plotted  on  a  semi-log  scale  the 

confidence  limits  will  be  centered  on  each  point. 

A  plot  of  f  as  a  function  of  degrees  of  freedom  k  for  several 

different  confidences  is  shown  in  Figure  9.  The  ordinate  expresses 

f  in  Db:  to  determine  f  one  goes  across  to  the  straight  line  labeled 
c  c 

"factor”,  then  reads  f  off  the  abscissa.  As  an  example,  for  k  =  30 


% 


the  90^  confidence  factor  expressed  in  Db  is  seen  to  be  «  3.8,  for 
which  f 

c 

The  next  two  sections  will  consider  formulas  for  determining 
k  for  each  of  the  two  methods  of  computing  PSDs. 

9.  k  for  the  MLP  Calculation 

The  number  of  degrees  of  freedom  k  for  the  one-dimensional  MLP 
calculated  PSD  is  given  by  (Blackman  and  Tukey,  1958) 


where 

N  =  the  number  of  data  points  in  the  data  set,  and 

M  =  the  maximum  number  of  lags  in  the  autocorrelation  function 
=  the  number  of  points  in  the  computed  PSD. 

As  can  be  seen  from  this  formula,  k  is  nearly  proportional  to 
the  ratio  N/M.  In  order  to  achieve  maximum  confidence  in  the  PSD 
estimates  it  is  necessary  that  k  be  made  as  large  as  possible;  this 
will  insure  that  the  confidence  factor  fc  [Eq.  (90)]  will  be  small. 

N  is  usually  some  fixed  number  of  points,  so  M  is  therefore  chosen  to 
be  small  relative  to  N.  It  is  for  this  reason  that  M  is  normally 
chosen  to  be  not  more  than  10—90^  of  N. 


k  for  the  FFT  Calculation 


7 

>• 


The  number  of  degrees  of  freedom  k  for  the  one-dimensional 
FFT  calculated  PSD  is  given  by  (Tukey,  19^7) 


i«  -  »,] 

N  ~  J 

o  ' 


where 

N  =  the  number  of  points  in  the  original  data  set, 

N,p  =  the  number  of  points  tapered  by  the  data  windov  g(t), 

Nq  =  the  total  number  of  data  points  after  zeros  have  been 
appended,  and 

p  =  the  effective  number  of  points  involved  in  smoothing  of 
the  raw  PSD  estimates. 

Although  the  MLP  and  FFT  formulas  for  k  appear  to  be  super¬ 
ficially  different,  it  is  not  difficult  to  show  that  the  confidence 
limits  calculated  from  a  given  k  are  essentially  equivalent  for  the 
two  cases,  as  would  be  expected. 

C.  Planning  Considerations 

In  the  previous  sections  the  basis  of  a  PSD  calculation  was 
discussed  with  appropriate  equations  given.  In  a  practical  applica¬ 
tion  some  thought  must  be  given  to  the  problem  of  choosing  data  sample 
spacing,  frequency  resolution,  and  the  frequency  range  over  which  the 
power  spectrum  is  to  be  calculated.  This  section  will  deal  only  with 
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considerations  to  be  made  in  performing  an  MLP  calculation;  similar 
considerations  will  apply  to  an  FFT  calculation. 

1.  Choosing  At  and  f 

°  max 

The  raw  data  from  which  the  PSD  is  to  be  calculated  is  assumed 

to  consist  of  discrete,  equi-spaeed  samples  of  period  At.  Then  if 

f  denotes  the  maximum  frequency  for  which  PSD  estimates  are  ob- 
max 

tained  we  have  from  the  sampling  theorem 


f  =  -i-  . 

max  PAt 

f  is  clearly  independent  of  the  total  number  of  data  points  and 
max 

the  maximum  lag  M  of  the  autocorrelation  function.  It  depends  solely 
on  the  sample  spacing.  If  in  an  analysis  one  wishes  to  examine  the 
PSD  up  through  a  specific  frequency,  this  formula  shows  what  the 
corresponding  minimum  sampling  frequency  must  be.  Likewise,  if  it  is 
desired  only  to  compute  the  PSD  for  a  limited  frequency  range,  by 
choosing  the  correct  At  the  resulting  PSD  will  have  a  frequency  range 
of  exactly  the  right  size.  This  fact  is  helpful  in  minimizing  the 
number  of  data  points  required  in  a  calculation.  For  example,  if  one 
were  interested  in  examining  the  PSD  of  a  set  of  data  only  in  the  fre¬ 
quency  range  0—5  Hz,  the  sampling  period  At  would  be  0.1  sec.  Samp¬ 
ling  more  often  than  this  would  increase  f  beyond  the  range  of 

max 

interest,  increase  the  total  number  of  data  points,  and  thus  the 
machine  calculation  time,  and  generally  would  not  add  any  information. 


If  the  sampled  data  already  exists  for  which  the  power  spectrum 
is  to  be  calculated,  it  is  possible  to  adjust  At  to  near  the  desired 
value  by  (a)  decimation,  using  only  every  p'th  point  of  the  original 
data  or  (b)  averaging,  averaging  by  groups  of  p  points.  The  former 
method  increases  At  by  a  factor  of  p  and  maintains  the  stendard  devia¬ 
tion  (for  random  or  new-random  data)  about  the  mean.  The  latter 
increases  At  by  a  factor  of  p  but  decreases  the  standard  deviation 
about  the  mean  by  a  factor  of  /p.  For  further  discussion  on  these 
methods  consult  Enochson  and  Otnes  (1968). 

2.  Choosing  M  and  Af 

The  freouency  resolution  Af  of  the  calculated  PSD  is 


where  M  is  the  maximum  lag  in  the  autocorrelation  fVinction.  This 
formula  may  be  used  for  either  determining  in  advance  of  a  calculation 
what  the  frequency  resolution  will  be,  or  for  determining  M  for  a 
given  desired  Af.  As  discussed  in  section  B,  M  should  be  small 
relative  to  the  total  amount  of  data.  This  fact  must  also  be  taken 
into  account  when  planning  a  PSD  calculation. 

As  an  example  of  how  the  foregoing  considerations  may  be  used, 
suppose  we  wish  to  calculate  the  PSD  of  a  data  set  for  the  frequency 


kO 


range  0—5  Hz.  Suppose  further  that  we  wish  the  frequency  resolution 
if  to  be  0.?  Hz,  providing  a  total  of  26  points  (including  the  end¬ 
points)  in  the  power  spectrum,  and  that  the  90^  confidence  factor  is 
to  be  <  2.0.  Then 

At  =  =0.1  sec  . 

M  ■  - ?-5  sec 

=  the  maximum  lag  of  26  data  points. 

From  Figure  9  we  find  that  for  the  desired  confidence  we  must 
have  k  >  45,  from  which 

t  >  f  .  „  fa  sec 

for  which 

N  *  6l0  data  points. 

As  a  second  example  suppose  a  data  set  already  exists  consist¬ 
ing  of  500  points  with  sample  spacing  1  sec.  If  the  autocorrelation 
function  is  truncated  at  lO’Jt  the  length  of  the  data  set,  we  calculate 

f _  .  Af ,  k,  and  f  (90^)  to  be: 

max  c 


4l 


f 

max 


1 

?At 


0.5  Hz 


9 


Af 


_1_  1 

2M  ~  ?(50) 


0.01  Hz 


k  = 


,  and 


fc  (90%)  = 2.9 


From  the  above  simple  examples  it  is  seen  that  M,  N,  At.  f  , 

Af,  and  f  are  to  a  degree  interrelated.  Therefore  one  must  be  some- 
max 

what  Judicious  in  their  specification  in  order  to  achieve  the  maximum 
amount  of  usable  information  in  a  PSD  calculation.  For  example,  one 
can  achieve  a  very  high  degree  of  statistical  reliability  by  making 
M  very  small.  However  this  would  be  at  the  expense  of  the  frequency 
resolution.  In  general,  for  a  given  data  set  there  is  a  trade-off 
between  statistical  reliability  and  frequency  resolution.  This  prob¬ 
lem  can  be  surmounted  by  simply  increasing  N  fusing  more  data)  but 
then  the  computing  time  increases,  and  so  on. 

D.  Averaging  Power  Spectra 

It  is  occasionally  desired  to  study  the  average  characteristics 
of  power  spectra  over  long  periods  of  time.  Although  it  is  possible 
to  compute  an  average  power  spectrum  over  a  long  data  set,  it  is 


sometimes  more  efficient  to  break  it  up  into  shorter  sets,  compute  the 
PSD  for  each  individual  data  set,  then  average  the  results.  This  pro¬ 
cedure  is  also  necessary  if,  for  example,  several  short,  non-contiguous 
data  sets  of  varying  lengths  exist  for  which  one  wishes  to  extract  an 
average  power  spectrum.  The  validity  of  the  averaged  results  will 
depend  on  the  stationarity  of  the  data;  for  data  which  is  non¬ 
stationary  the  results  will  be  influenced  by  the  lengths  and  number  of 
data  sets  used  in  the  calculation.  [For  a  short  discussion  on  this 
point  see  Sentman  (19?^).  1  The  following  sections  apply  to  MLP  cal¬ 
culated  PSD's,  though  similar  considerations  would  apply  to  FFT  cal¬ 
culated  PSD's. 

1.  Averaged  Pbwer  Spectra  from  Data  Sets  of  Equal  Lengths 

If  it  is  assumed  that  the  frequency  resolution  Af  is  identical 
for  each  PSD  comprising  the  average,  the  statistical  reliability  of 
each  is  the  same.  Then 

?i  *  5  £  PU 

where 

Pj  =  i'th  point  in  the  averaged  power  spectrum, 

=  i'th  point  in  the  j'th  power  spectrum, 

N  =  the  number  of  spectra  being  averaged,  and 


h'j 


<3j  =  standard  deviation  about  the  mean  of  data  in  the  original 

j'th  (pre- whitened)  data  set.  If  the  data  in  the  j'th 
data  set  were  divided  by  a  after  the  mean  was  removed 

J 

prior  to  calculation  of  the  power  spectrum,  this  weight 
is  necessary  to  restore  the  proper  relative  units  to 
If  the  data  were  not  normalized  this  weight  is 
equal  to  one. 

?.  Averaged  Power  Spectra  from  Data  Sets  of  Different  Lengths 
If  the  individual  spectra  comprising  the  average  are  computed 
from  data  sets  of  different  lengths,  the  statistical  reliability  of 
each  PSD  is  different.  If  Af  is  identical  for  each  individual  PSD. 
the  expression  for  averaging  becomes 


"  -  b  I 
?Q/kj  -  lj 


where 

k.  =  degrees  of  freedom  in  the  j'th  computed  power  spectrum. 

J 

and 

b  =8,  for  50^  confidence  limits. 

The  exponential  weighting  factor  must  be  included  to  properly 
weigh  individual  spectra  according  to  its  probable  error.  The  uncer¬ 
tainty  in  the  spectra  is  taken  to  be  the  probable  error,  or  the 
square  root  of  the  50j6  confidence  factor  defined  in  Eq.  (20). 
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Statistical  Reliability  of  Averaged  Spectra 
Confidence  limits  for  the  averaged  spectrum  are  computed  by 
assuming  that  its  equivalent  degrees  of  freedom  it  equal  the  sum  of 
the  degrees  of  freedom  in  the  individual  spectra  comprising  the  aver¬ 
age. 


k  = 


It  can  easily  be  shown  that  if  the  statistical  "noise"  present  in  each 
of  the  individual  spectra  is  treated  as  a  random  fluctuation  about 
the  true  power  spectrum,  the  above  expression  for  calculating  k  results 
in  confidence  limits  that  shrink  at  exactly  the  same  rate  as  the 
"noise"  when  more  and  more  spectra  are  averaged. 


V.  AN  EXAMPLE  OF  A  POWER  SPECTRUM  CALCULATION 


As  an  example  of  how  the  information  contained  in  the  previous 
sections  may  be  applied,  the  following  will  serve  to  illustrate  the 
basic  features  of  an  MLP  calculated  power  spectrum.  The  example 
given  is  from  a  calculation  made  in  conjunction  with  the  study  of 
oscillatory  phenomena  in  the  solar  atmosphere  (Sentman,  197'). 

The  raw  data  consisted  of  antenna  temperature  5  sec  data 
samples  of  solar  microwave  emission  recorded  on  the  North  Liberty 
Radio  Observatory  2-cm  radiometer.  Individual  data  sets  ranged  in 
length  from  k — 12  h.  The  power  spectrum  of  each  data  set  was  to  be 
calculated,  and  the  average  of  all  the  individual  spectra  computed 
to  obtain  an  average  power  spectrum  for  all  the  data  sets. 

The  frequency  range  of  interest  to  this  study  was  0—15  mHz 

_  T 

(l  mHz  =  10- '  Hz),  or  f  =15  mHz.  A  sampling  tinie  At  of 

max 

l/2(l5  x  10  )  =  53,5  sec  was  thus  indicated.  The  nearest  that  one 

could  come  to  this  figure  using  5  sec  data  samples  was  either  ^0  sec 

or  35  aec,  corresponding  to  averaging  the  5  sec  samples  by  groups  of 

6  or  7,  respectively.  Averaging  by  groups  of  6  was  chosen,  yielding 

an  effective  sampling  time  At  =  50  sec  an  i  a  maximum  frequency 

f  =16.7  mHz.  A  typical  data  set  with  30  sec  resolution  (effective 
max 

sampling  time)  is  shown  in  Figure  10. 


It  was  known  in  advance  of  the  calculation  (by  trial  runs  on 
several  data  sets)  that  the  power  spectra  all  contained  a  very  large 
peak  at  near  zero  frequencies,  thus  indicating  a  need  for 
pre-whitening.  High-pass  smoothing  was  achieved  by  means  of  the 
method  described  in  section  III-B  by  making  10  passes  through  the 
data  with  the  smoothing  function  described  by  Eq.  (12).  The  result¬ 
ing  low-  and  high-pass  power  transfer  functions  are  shown  in  Figure 
11.  A  maximum  lag  of  25  min  (=  50  lags  *  50  sec)  was  chosen  to 
balance  frequency  resolution  against  statistical  reliability.  This 
resulted  in  a  frequency  resolution  Af  =  0.55  mHz.  The  degrees  of 
freedom  k  therefore  ranged  from  18.5  (4  h  data  set)  to  5^.9  (12  h 
data  set),  providing  the  high  degree  of  statistical  reliability 
necessary  for  the  study  (very  low  amplitude  fluctuations  were  being 
sought ) . 

The  data  normalization  and  power  spectra  calculations  were 
achieved  using  the  subroutine  listed  in  Appendix  I.  After  the  power 
spectrum  was  calculated  it  was  post-darkened  to  compensate  for  the 
pre-whitening  by  multiplying  with  the  inverse  of  the  high-pass  power 
transfer  function  shown  in  Figure  11.  An  example  of  the  resulting 
spectra  with  confidence  limits,  normalized  to  a  value  of  10  in  the 
range  of  4— 5  mHz  for  display  purposes,  ds  shown  in  Figure  12  for 
the  frequency  range  0—15  mHz.  All  calculations  were  carried  out  on 
a  Univac  4l8  computer. 


FIGURE  CAPTIONS 


Figure  9  The  confidence  factor  as  a  function  of  degrees  of  freedom 
k  (see  page  35). 

Figure  10  Typical  plot  of  antenna  temperature  versus  time  with  30 
sec  time  resolution.  Data  records  used  to  compute  power 
spectra  were  chosen  to  exclude  the  end  pieces  of  each  data 


day. 

Figure  11  Low-pass  and  high-pass  filter  functions  R(f)  and  R*(f) 
used  to  pre-whiten  the  data  prior  to  calculation  of  the 
power  spectra.  Post -darkening  is  achieved  by  multiplying 
the  resultant  spectra  by  l/R#(f). 

Figure  12  Power  spectrum  of  the  single  record  containing  a  statis¬ 
tically  significant  peak  near  ^  mHz. 
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APPENDIX  I 

A  FORTRAN  SUBROUTINE  FOR  CALCULATING  POWER  SPECTRA 

The  following  is  an  example  of  a  subroutine  for  an  MLP  calcu¬ 
lation  of  a  power  spectrum.  Pre-whitening  is  assumed  to  have  been 
coepleted  prior  to  entry  into  the  subroutine,  and  post-darkening  and 
conversion  to  absolute  units  from  normalized  units  is  assumed  to 
take  place  after  return. 

The  program  was  written  from  the  definitions  for  the  auto¬ 
correlation  function  and  discrete  cosine  transform.  However,  to 
facilitate  the  transform  operation  a  cosine  table  is  constructed  and 
a  table  look-up  procedure  is  incorporated  into  the  program  rather 
than  requiring  a  double  precision  cosine  to  be  computed  each  time 
it  is  needed.  The  table  is  constructed  in  the  first  call  of  the 
subroutine  and  each  time  the  autocorrelation  function  maximum  lag  is 
changed;  succeeding  calls  use  the  table  already  constructed. 

As  shown  below,  the  program  will  accommodate  a  maximum  of  1500 
data  points  with  an  autocorrelation  function  maximum  lag  of  500. 

These  figures  may  be  raised  or  lowered  as  desired  by  dimensioning  the 
relevant  arrays  accordingly. 
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A.  Subroutine  Arguments 

Input  variables: 

X  -  Double  precision  array  containing  data  to  be  power 
spectrum  analyzed 

N  -  Integer  specifying  the  number  of  data  points  In  X  (may 
be  smaller  than  the  dimension  size  of  X) 

MA  -  Integer  specifying  the  maximum  number  of  lags  for  the 

autocorrelation  function  (may  be  smaller  than  the  dimen¬ 
sion  size  of  R) 

Output  variables: 

U  -  Single  precision  array  containing  M  +  1  hamming  smoothed 
normalized  power  spectral  estimates 

R  -  Double  precision  array  containing  M  +  1  values  of  the 
normalized  autocorrelation  coefficients 

XBAR  -  Double  precision  variable  for  the  mean  of  the  input 
data  set 

SX  -  Double  precision  variable  for  the  standard  deviation 
of  the  input  data  set.  This  variable  is  necessary  if 
absolute  units  are  to  be  calculated  for  the  power  spec¬ 
trum  (not  done  here) 
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B.  Subroutine 


SUBROUTINE  SPBCTR(X,N,MA,U,R,XBAR,SX) 

DOUBLE  PRECISION  X,AL,CS,XBAR,SX,R 

DOUBLE  PRECISION  FAC,ARG,DCOS,DBLE,DSQRT,PI,AEND 

DIMENSION  X  (1500 )  ,U(501 )  ,R  (501 ) ,  CS  (501 ) ,  AL  (501 ) 

DATA  PI/.3i415O2654D+01/,M/0/ 

C 

C  NORMALIZE  DATA  TO  ZERO  MEAN  AND  UNIT  STANDARD  DEVIATION 

XBAR-0.0 
DO  3  1*1. N 
3  XBAR-XBAR+XU) 

XBAR-XBAR/DBLE  (FLOAT  (N  ) ) 

SX-O.OD+OO 
DO  6  1*1. N 
X(I)-X(I)-XBAR 

6  SX«SX+X(l)*X(l) 

SXaDSQRT  (sx/dble  (float  (n  ) ) ) 

DO  7  I-l.N 

7  X(I)«X(I)/SX 
C 

C  BUILD  COSINE  TABLE  .CS.  FOR  INTERVAL  .0.  to  .PI. 

IF(MA.EQ.M)  00  TO  200 
M44A 
MP^+1 

AVERB-DBLE  (FLOAT  (M  ) ) 

DO  2  I-1,MP 
ARO-DBLE  (FLOAT  (i-l )  ) 

ARC*  (  PI*ARG  )  /AVERS 
2  CS ( I )*DCOS  (ARO ) 

200  CONTINUE 
C 

C  CALCULATE  AUTOCORRELATION  COEFFICIENTS 

DO  30  J-1,MP 
r(j )*o.od+oo 

JP-J-1 

IEND-N-JP 

AEND*DBI£  (FLOAT  (  TEND  ) ) 

DO  25  I-1,IEND 
INDX-I+JP 

25  R^i-RUj+Xd^XflNPX) 

30  R(J)-R(J)/AHn> 

C 


ooo  nan  o  n 
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CALCULATE  POWER  SPECTRUM 
(FOURIER  TRANSFORM  A/C  FUNCTION  .rt. ) 

MX?W2 

MBC240C2+2 

MM4-1 

DO  40  J«1,MP 

AL(j)-0.0D+OO 

JP-J-1 

DO  55  K-1,MM 

CALCULATE  INDEX  .INDCS.  FOR  RETRIEVING 
PROPER  COSINE  FROM  TABLE 
INDCS^OD(  (JP*K )  ,MX2  )+l 
IF ( INDCS. GT.MP)  INDCSaMPX2- INDCS 
55  AL(J)«AL(J)+2.0D+00*R(K+1)*CS (INDCS) 

FAC-1. OD+OO 

IF(M0D(JP,2).EQ.l)  FAC— l.OD+OO 
UO  AL  ( J  )-AL  ( J  )+R  ( 1 )  +R  (MP  )*FAC 

APPLY  HAM4ING  SMOOTHING  FUNCTION  TO  ESTIMATES  .AL. 
TO  YIELD  SMOOTHED  POWER  SreCTRUM  ESTIMATES  .U. 

U(l  )-SNGL  (0. 5UD+00*AL  (l  )+0.  U6D+00»Al(2  )  ) 

U  (MP)-SNGL  (0. 5UD+00»AL  (MP  )+0.  U6D+00*AL  (M  )  ) 

DO  50  I«2,M 

50  U(I)-SNBL(0.51*D+00»AL(I)+0.25D^X)»(AL(I-1)+AL(I+1))) 
RETURN 
END 
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